
#include<stdio.h>
#include<stdlib.h>
#include<mpi.h>
#include "comm.h"
#include "allocate.h"
#include "rtmlib.h"
#include "gptl.h"

int main(int argc,char**argv){
  //初始化mpi环境
  MPI_Init(&argc,&argv);
  ret=GPTLinitialize();
  gs("Total");
  gs("Prepation");
  //线程库初始化
  athread_init(); 
  //获取进程数及进程id
  MPI_Comm_size(MPI_COMM_WORLD,&process_size);
  MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
  //进程参数
  MPI_Comm cartcomm;
  dims[0]=atoi(argv[1]);
  dims[1]=atoi(argv[2]);
  dims[2]=atoi(argv[3]);
  int periods[3]={0,0,0},reorder=0,coords[3];
  //创建拓扑进程
  MPI_Cart_create(MPI_COMM_WORLD,3,dims,periods,reorder,&cartcomm);
  //进程行\列坐标
  MPI_Cart_coords(cartcomm,my_rank,3,coords);
  x_rank=coords[0];
  y_rank=coords[1];
  z_rank=coords[2];
  //获得当前进程的邻居 
  MPI_Cart_shift(cartcomm,0,1,&nbrs[X0],&nbrs[X1]);
  MPI_Cart_shift(cartcomm,1,1,&nbrs[Y0],&nbrs[Y1]); 
  MPI_Cart_shift(cartcomm,2,1,&nbrs[Z0],&nbrs[Z1]); 
  Init_parmter();
  //模型初始化及扩展模型-边界吸收
  if(my_rank==0){
	Init_PML();
	vpmodel=allocate1float((mgridnumx*mgridnumy*mgridnumz));
	vp_ext=allocate1float(((mgridnumx+2*PMLX)*(mgridnumy+2*PMLY)*(mgridnumz+2*PMLZ)));
	Init_model();
	generate_v_3d(vpmodel,mgridnumx,mgridnumy,mgridnumz, PMLX,PMLY,PMLZ, vp_ext);	
	/*check extend model */
	//FILE*fp =fopen("./model/model_ext.dat","wb");
	//fwrite(vp_ext,sizeof(float),((mgridnumx+2*PMLX)*(mgridnumy+2*PMLY)*(mgridnumz+2*PMLZ)),fp);
	////fwrite(vpmodel,sizeof(float),((mgridnumx)*(mgridnumy)*(mgridnumz)),fp);
	//fclose(fp);
	//
  }
  //更新的dt,Vpmax,PMLX,PMLZ分发所有进程
  MPI_Bcast(&dt,1,MPI_FLOAT,0,MPI_COMM_WORLD);
  MPI_Bcast(&Vpmax,1,MPI_FLOAT,0,MPI_COMM_WORLD);
  MPI_Bcast(&Vpconst,1,MPI_FLOAT,0,MPI_COMM_WORLD);
  MPI_Bcast(&PMLX,1,MPI_INT,0,MPI_COMM_WORLD);
  MPI_Bcast(&PMLY,1,MPI_INT,0,MPI_COMM_WORLD);
  MPI_Bcast(&PMLZ,1,MPI_INT,0,MPI_COMM_WORLD);
  //任务划分
  if((mgridnumx+2*PMLX)%dims[0]==0&&(mgridnumy+2*PMLY)%dims[1]==0&&(mgridnumz+2*PMLZ)%dims[2]==0){
	//提示任务划分的信息
	DivNx=(mgridnumx+2*PMLX)/dims[0];
	DivNy=(mgridnumy+2*PMLY)/dims[1];
	DivNz=(mgridnumz+2*PMLZ)/dims[2];
	SubNx=DivNx+10;//通信区域为5
	SubNy=DivNy+10;
	SubNz=DivNz+10;
	if(my_rank==0){
	  if(DivNx<PMLX){
		printf("X方向:子区域DivNx=%d小于吸收层厚度=%d，请减少该方向进程数!\n",DivNx,PMLX);
		fprintf(flog,"X方向:子区域DivNx=%d小于吸收层厚度=%d，请减少该方向进程数!\n",DivNx,PMLX);
	  }
	  if(DivNy<PMLY){
		printf("Y方向:子区域DivNy=%d小于吸收层厚度=%d，请减少该方向进程数!\n",DivNy,PMLY);
		fprintf(flog,"Y方向:子区域DivNy=%d小于吸收层厚度=%d，请减少该方向进程数!\n",DivNy,PMLY);
	  }
	  if(DivNz<(PMLZ+rz)){
		printf("Z方向:子区域DivNz=%d小于吸收层厚度=%d，请减少该方向进程数!\n",DivNz,PMLZ+rz);
		fprintf(flog,"Z方向:子区域DivNz=%d小于吸收层厚度=%d，请减少该方向进程数!\n",DivNz,PMLZ+rz);
	  }
	}
  }
  else{
	if(my_rank==0){
	  printf("当前进程数无法整除!\n");
	  fprintf(flog,"当前进程数无法整除!\n");
	}
	MPI_Finalize();
	return 0;
  }
  if(my_rank == 0){
	printf("Max Vp velocity is %.3fm/s, Suggested time interval dt=%.6f\n",Vpmax,dt);
	printf("Original model size      :mgridnumx=%d,mgridnumy=%d,mgridnumz=%d\n",mgridnumx,mgridnumy,mgridnumz);
	printf("Suggested absorbed layer :PMLX=%d,PMLY=%d,PMLZ=%d\n",PMLX,PMLY,PMLZ);
	printf("Total calculation size   :NX=%d,NY=%d,NZ=%d\n",mgridnumx+2*PMLX,mgridnumy+2*PMLY,mgridnumz+2*PMLZ);
	printf("SubDomain size           :SubNx=%d,SubNy=%d,SubNz=%d \n",SubNx,SubNy,SubNz);
	printf("Effect size              :DivNx=%d,DivNy=%d,DivNz=%d \n",DivNx,DivNy,DivNz);

	fprintf(flog,"Max Vp velocity is %.3f, Suggested time interval dt=%.6f\n",Vpmax,dt);
	fprintf(flog,"Original model size      :mgridnumx=%d,mgridnumy=%d,mgridnumz=%d\n",mgridnumx,mgridnumy,mgridnumz);
	fprintf(flog,"Suggested absorbed layer :PMLX=%d,PMLY=%d,PMLZ=%d\n",PMLX,PMLY,PMLZ);
	fprintf(flog,"Total calculation size   :NX=%d,NY=%d,NZ=%d\n",mgridnumx+2*PMLX,mgridnumy+2*PMLY,mgridnumz+2*PMLZ);
	fprintf(flog,"SubDomain size           :SubNx=%d,SubNy=%d,SubNz=%d \n",SubNx,SubNy,SubNz);
	fprintf(flog,"Effect size              :DivNx=%d,DivNy=%d,DivNz=%d \n",DivNx,DivNy,DivNz);
	fflush(flog);
  }
  //分发模型
  distribute_model();
  //生成震源子波
  wave_len=generate_wavelet();
  //炮循环
  //申请波场数据
  allocate_wavefield();
  //申请通信buffer
  allocate_buffer();
  //初始化震源位置
  int* SX=allocate1int(Nsx);
  int* SY=allocate1int(Nsy);
  int  SZ=sz;
  int ix,iy;
  for(ix=0;ix<Nsx;ix++){
	SX[ix]=(ix+1)*((mgridnumx-1)*1.0f/(Nsx+1));
  }
  for(iy=0;iy<Nsy;iy++){
	SY[iy]=(iy+1)*((mgridnumy-1)*1.0f/(Nsy+1));
  }
  if(my_rank==0){
	printf("Total shot X,Y is %d,%d \n",Nsx,Nsy);
	fprintf(flog,"Total shot X,Y is %d,%d \n",Nsx,Nsy);
	fflush(flog);
  }
  double t1,t2,t3,t4,t5,t6,t7,t8;//计时器
  Init_kernel_par();
  //初始化总能量
  memset(Total_image,0,sizeof(float)*SubNx*SubNy*SubNz);
  memset(shot_energy,0,sizeof(float)*SubNx*SubNy*SubNz);
  ge("Prepation");
  gs("Total Migration");
  t7=MPI_Wtime();
  for(ix=0;ix<Nsx;ix++){
	for(iy=0;iy<Nsy;iy++){
	  gs("Each Migration");
	  gs("Each First Arrival modeling");
	  t5=MPI_Wtime();
	  sx=SX[ix]+PMLX;
	  sy=SY[iy]+PMLY;
	  sz=SZ+PMLZ;
	  if(my_rank==0)  printf("%d ,%d \n",sx,sy);
	  
	  if(my_rank==0)  printf("memset 0 start\n");
	  gs("memset");
	  //内存置零
	  memset(T,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Txx,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Tyy,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Tzz,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Vx ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Vy ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Vz ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(image ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(recordVz,0, (DivNx*DivNy*Nt)*sizeof(float));
	  ge("memset");
	  if(my_rank==0)  printf("memset 0 end\n");
	  //模拟直达波切除
	  parameter.Vpconst=Vpconst;//设置直达波模式
	  for(it=0;it<Nt;it++){
		if(my_rank==0&&it%50==0){
		  fprintf(flog,"Current shot is x-y: %d-%d ,shot position grid x,y,z : %d,%d,%d \n",ix,iy,sx,sy,sz);
		  fprintf(flog,"First Arrival modeling step is %d \n",it);
		  printf("Current shot is x-y: %d-%d ,shot position grid x,y,z : %d,%d,%d \n",ix,iy,sx,sy,sz);
		  printf("First Arrival modeling is %d \n",it);
		  fflush(stdout);
		}
		gs("first");
		/*1. merge sress t= txx+tyy+tzz*/
		gs("mergeStress-fr");
		mergeStress();
		ge("mergeStress-fr");
		t1=MPI_Wtime();
		/*2.update velocity  */
		gs("Cal-V1-fr");
		update_Vel1();
		ge("Cal-V1-fr");

		gs("mpi-V-fr");
		MPI_LaunchV();
		ge("mpi-V-fr");

		gs("Cal-V2-fr");
		update_Vel2();
		athread_join(); //回收从核
		ge("Cal-V2-fr");


		gs("mpiV-join-fr");
		MPI_JoinV();
		ge("mpiV-join-fr");
		if((x_rank==sx/DivNx)&&(y_rank==sy/DivNy)&&(z_rank==sz/DivNz)){//确定震源位置所在进程
		  add_souce((sx-(sx/DivNx)*DivNx+5),(sy-(sy/DivNy)*DivNy+5),(sz-(sz/DivNz)*DivNz+5),it);//重新更正震源坐标
		}
		t2=MPI_Wtime();

		gs("save_record-fr");
		save_record(it);  
		ge("save_record-fr");

		t3=MPI_Wtime();
		gs("Cal-T1-fr");
		/*3.update stress */
		update_T1();//更新应力场
		ge("Cal-T1-fr");

		MPI_LaunchT();
		gs("Cal-T2-fr");
		update_T2();
		athread_join();//回收从核
		ge("Cal-T2-fr");

		gs("mpi-T-join-fr");
		MPI_JoinT();
		ge("mpi-T-join-fr");
		t4=MPI_Wtime();
		ge("first");
		if(my_rank==0&&it%50==0){
		  printf("Each step: Vel=%lf,T=%lf\n",(t2-t1),(t4-t3));
		  fprintf(flog,"Each step: Vel=%lf,T=%lf\n",(t2-t1),(t4-t3));
		  fflush(flog);
		  fflush(stdout);
		}
		//outforwardsnap(it,ix,iy);//输出快照
	  }
	  t6=MPI_Wtime();
	  if(my_rank==0)  printf("Arrival  modeling runtime=%lf\n",(t6-t5));
	  if(my_rank==0)  fprintf(flog,"Arrival  modeling runtime=%lf\n",(t6-t5));
	  ge("Each First Arrival modeling");
#if 0
	  gs("Each shot modeling");
	  t5=MPI_Wtime();
	  //正演
	  memset(T,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Txx,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Tyy,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Tzz,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Vx ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Vy ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(Vz ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  parameter.Vpconst=0;//设置正常模拟
	  for(it=0;it<Nt;it++){
		if(my_rank==0&&it%50==0){
		  fprintf(flog,"Current shot is x-y: %d-%d ,shot position grid x,y,z : %d,%d,%d \n",ix,iy,sx,sy,sz);
		  fprintf(flog,"Forwad  step is %d \n",it);
		  printf("Current shot is x-y: %d-%d ,shot position grid x,y,z : %d,%d,%d \n",ix,iy,sx,sy,sz);
		  printf("Forwad  step is %d \n",it);
		  fflush(flog);
		  fflush(stdout);
		}
		gs("each shot");
		/*1. merge sress t= txx+tyy+tzz*/
		gs("merge");
		mergeStress();
		ge("merge");

		t1=MPI_Wtime();
		/*2.update velocity  */
		gs("Cal-V");
		gs("Cal-V1");
		update_Vel1();
		ge("Cal-V1");

		gs("boundary-Pack-T");
		boundary_packT(it);	
		ge("boundary-Pack-T");

		gs("mpi-Pack-V");
		MPI_LaunchV();
		ge("mpi-Pack-V");

		gs("Cal-V2-launch");
		update_Vel2();
		ge("Cal-V2-launch");

		gs("boundary-write-T");
		if((it+1)%Nf==0)
		  boundary_writeT(it);
		ge("boundary-write-T");

		gs("Cal-V2-join");
		athread_join(); //回收从核
		ge("Cal-V2-join");

		gs("mpi-V-join");
		MPI_JoinV();
		ge("mpi-V-join");

		if((x_rank==sx/DivNx)&&(y_rank==sy/DivNy)&&(z_rank==sz/DivNz)){//确定震源位置所在进程
		  add_souce((sx-(sx/DivNx)*DivNx+5),(sy-(sy/DivNy)*DivNy+5),(sz-(sz/DivNz)*DivNz+5),it);//重新更正震源坐标
		}
		ge("Cal-V");
		t2=MPI_Wtime();


		gs("Cal-T");

		gs("save_record");
		save_record(it);  
		ge("save_record");

		t3=MPI_Wtime();
		/*3.update stress */

		gs("Cal-T1");
		update_T1();//更新应力场
		ge("Cal-T1");

		gs("boundary_packV");
		boundary_packV(it);	
		ge("boundary_packV");

		gs("mpi-Pack-T");
		MPI_LaunchT();
		ge("mpi-Pack-T");

		gs("Cal-T2-launch");
		update_T2();
		ge("Cal-T2-launch");


		gs("boundary_write-V");
		if((it+1)%Nf==0)
		  boundary_writeV(it);
		ge("boundary_write-V");

		gs("Cal-T2-join");
		athread_join();//回收从核
		ge("Cal-T2-join");


		gs("mpi-T-join");
		MPI_JoinT();
		ge("mpi-T-join");

		ge("Cal-T");
		t4=MPI_Wtime();
		if(my_rank==0&&it%50==0){
		  printf("Each step :Vel=%lf,T=%lf\n",(t2-t1),(t4-t3));
		  fprintf(flog,"Each step :Vel=%lf,T=%lf\n",(t2-t1),(t4-t3));
		  fflush(flog);
		  fflush(stdout);
		}

		/*4. output snap */
		gs("outforwardsnap");
		outforwardsnap(it,ix,iy);//输出快照
		ge("outforwardsnap");

		ge("each shot");
	  }
	  //保存当前炮记录 
	  //write_record(ix,iy);
	  t6=MPI_Wtime();
	  if(my_rank==0)  printf("shot  modeling runtime=%lf\n",(t6-t5));
	  if(my_rank==0)  fprintf(flog,"shot  modeling runtime=%lf\n",(t6-t5));
	  ge("Each shot modeling");
	  gs("Each backward & image");
	  t5=MPI_Wtime();
	  //逆时成像
	  memset(extrapT,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(extrapTxx,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(extrapTyy,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(extrapTzz,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(extrapVx ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(extrapVy ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  memset(extrapVz ,0,(SubNx*SubNy*SubNz*sizeof(float)));	
	  for(it=Nt-1;it>wave_len/2;it--){
		if(my_rank==0&&it%50==0){
		  printf("Current shot is x-y %d-%d ,shot position grid x,y,z: %d,%d,%d \n",ix,iy,sx,sy,sz);
		  printf("Backward  step is %d \n",it);                           
		  fprintf(flog,"Current shot is x-y %d-%d ,shot position grid x,y,z: %d,%d,%d \n",ix,iy,sx,sy,sz);
		  fprintf(flog,"Backward  step is %d \n",it);                     
		  fflush(flog);
		  fflush(stdout);
		}
		t1=MPI_Wtime();
		//成像处理
		gs("image");
		imaging();
		ge("image");
		t2=MPI_Wtime();

		t3=MPI_Wtime();
		gs("outbackward");
		outbackward_extrap(it,ix,iy);
		//输出当前成像结果
		outimage(it,ix,iy);
		ge("outbackward");
		t4=MPI_Wtime();
		if(my_rank==0&&it%50==0){
		  printf("Each step: Image =%lf, outbackward =%lf\n",t2-t1,t4-t3);
		  fprintf(flog,"Each step: Image =%lf, outbackward =%lf\n",t2-t1,t4-t3);
		  fflush(flog);
		  fflush(stdout);
		}
		/*波场重建*/
		/*1.*/
		t1=MPI_Wtime();

		gs("boundary-read-T");
		if((it+1)%Nf==0)
		  boundary_readV(it);
		ge("boundary-read-T");


		gs("boundary-unPack-V");
		boundary_unpackV(it);
		ge("boundary-unPack-V");

		gs("Cal-T1-back");
		update_T1_back();
		ge("Cal-T1-back");

		gs("mpi-Pack-T-back");
		MPI_LaunchT();
		ge("mpi-Pack-T-back");


		gs("Cal-T2-back-launch");
		update_T2_back();
		ge("Cal-T2-back-launch");


		gs("Cal-T2-back-join");
		athread_join();
		ge("Cal-T2-back-join");


		gs("mpi-T-back-join");
		MPI_JoinT();
		ge("mpi-T-back-join");

		t2=MPI_Wtime();
		t3=MPI_Wtime();
		/*2.*/

		gs("boundary-read-V");
		if((it+1)%Nf==0)
		  boundary_readT(it);
		ge("boundary-read-V");


		gs("boundary-unPack-T");
		boundary_unpackT(it);
		ge("boundary-unPack-T");


		gs("mergeStress");
		mergeStress();
		ge("mergeStress");

		gs("Cal-V1-back");
		update_Vel1_back();
		ge("Cal-V1-back");


		gs("mpi-Pack-V-back");
		MPI_LaunchV();
		ge("mpi-Pack-V-back");


		gs("Cal-V2-back-launch");
		update_Vel2_back();
		ge("Cal-V2-back-launch");

		gs("Cal-V2-back-join");
		athread_join();
		ge("Cal-V2-back-join");

		gs("mpi-V-back-join");
		MPI_JoinV();
		ge("mpi-V-back-join");

		t4=MPI_Wtime();
		if(my_rank==0&&it%50==0){
		  fprintf(flog,"Each step :T =%lf, Vel=%lf\n",t2-t1,t4-t3);
		  fflush(flog);
		  printf("Each step :T =%lf, Vel=%lf\n",t2-t1,t4-t3);
		  fflush(stdout);
		}

		/*波场外推*/
		t1=MPI_Wtime();

		load_record(it);

		gs("Cal-T1-extrap");
		update_T1_extrap();
		ge("Cal-T1-extrap");

		gs("mpi-T-extrap-Launch");
		MPI_LaunchT_extrap();
		ge("mpi-T-extrap-Launch");

		gs("Cal-T2-extrap-launch");
		update_T2_extrap();
		ge("Cal-T2-extrap-launch");

		gs("Cal-T2-extrap-join");
		athread_join();
		ge("Cal-T2-extrap-join");

		gs("mpi-T-extrap-join");
		MPI_JoinT_extrap();
		ge("mpi-T-extrap-join");

		t2=MPI_Wtime();
		t3=MPI_Wtime();

		gs("mergeStress_extrap");
		mergeStress_extrap();
		ge("mergeStress_extrap");

		gs("Cal-V1-extrap");
		update_Vel1_extrap();
		ge("Cal-V1-extrap");

		gs("mpi-V-extrap-Launch");
		MPI_LaunchV_extrap();
		ge("mpi-V-extrap-Launch");


		gs("Cal-V2-extrap-Launch");
		update_Vel2_extrap();
		ge("Cal-V2-extrap-Launch");


		gs("Cal-V2-extrap-join");
		athread_join();
		ge("Cal-V2-extrap-join");

		gs("mpi-V-extrap-join");
		MPI_JoinV_extrap();
		ge("mpi-V-extrap-join");
		t4=MPI_Wtime();
		if(my_rank==0&&it%50==0){
		  fprintf(flog,"Each step :T =%lf, Vel=%lf\n",t2-t1,t4-t3);
		  printf("Each step :T =%lf, Vel=%lf\n",t2-t1,t4-t3);
		  fflush(flog);
		  fflush(stdout);
		}
	  }
	  if(my_rank==0){
		printf("Stack current shot image. \n");
		fprintf(flog,"Stack current shot image. \n");
		fflush(flog);
		fflush(stdout);
	  }

	  gs("stack_imaging");
	  stack_image();
	  ge("stack_imaging");
	  if(my_rank==0){
		printf("Process total image. \n");
		fprintf(flog,"Process total image. \n");
		fflush(flog);
	  }

	  gs("post_imaging");
	  //成像后处理
	  if(ix==(Nsx-1)&&iy==(Nsy-1))
		post_imaging(ix,iy);	
	  ge("post_imaging");
	  //成像后处理
	  if(my_rank==0){
		printf("Current shot Finished. \n");
		fprintf(flog,"Current shot Finished. \n");
		fflush(flog);
	  }
	  t6=MPI_Wtime();
	  if(my_rank==0)  printf("each migration runtime=%lf\n",(t6-t5));
	  if(my_rank==0)  fprintf(flog,"each migration runtime=%lf\n",(t6-t5));
#endif
	  ge("Each backward & image");
	  ge("Each Migration");
	}
  }
  ge("Total Migration");
  if(my_rank==0){
	printf("All shot Finished. \n");
	fprintf(flog,"All shot Finished. \n");
	fflush(flog);
  }
  gs("Free memory");
  //释放空间
  free_buffer();
  free_wavefield();
  ge("Free memory");
  if(my_rank==0){
	printf("Free memomry Finished. \n");
	fprintf(flog,"Free memomry Finished. \n");
	fflush(flog);
  }
  t8=MPI_Wtime();
  if(my_rank==0)  printf("total  runtime=%lf\n",(t8-t7));
  if(my_rank==0)  fprintf(flog,"total runtime=%lf\n",(t8-t7));
  ge("Total");
	if(my_rank==0){
		sprintf(gptlfile,"./log/group%d.timingrank=%d",nsgroup,my_rank);
		GPTLpr_file(gptlfile);
	}
	//结束并行环境
  MPI_Finalize();
  return 0;
}

